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p I ■ Abstract 

■ Cells grown in culture act as a model system for analyzing the effects of anticancer compounds, 

. which may affect cell behavior in a cell cycle position-dependent manner. Cell synchronization 

techniques have been generally employed to minimize the variation in cell cycle position. However, 
synchronization techniques are cumbersome and imprecise and the agents used to synchronize the 
cells potentially have other unknown effects on the cells. An alternative approach is to determine 
the age structure in the population and account for the cell cycle positional effects post hoc. Here 
we provide a formalism to use quantifiable age distributions from live cell microscopy experiments 
^ | to parameterize an age-structured model of cell population response. 
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1 Introduction 

When examined individually in time lapse microscopy experiments, cells grown in culture display 
variability in the length of their cell cycles (between mitotic events), and this variability is represented 
by intermitotic time (IMT) distributions [201 l32l [33] . These distributions are usually obtained from 
asynchronously dividing populations of cells, which achieve a steady-state age structure when the 
population is growing exponentially. In experimental studies that examine the effects of perturbations 
on cellular proliferation, it is often desirable to determine whether the perturbation is affecting cells 
in particular stages within their cell cycle, i.e. whether the perturbation has cell cycle-specific effects. 
In this paper we provide a formalism to convert data obtained as IMT distributions to parameterize 
an age-structured population model and, thus, identifying the contribution of age structure to the 
response of the cell population to perturbation. 

Models of age-structured populations using partial differential equations, such as those originally 
discribed by A. Lotka, A.G. McKendrick, W.O. Kermack, and F. von Foerster, are well adapted to 
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model the dynamical features of experimental cultures of cells transiting the cell cycle with variable 
IMT. These models have been widely studied from a mathematical perspective [4"[ l3[ [T4 t [Pol \TE[ l2"T | l34l 
I3"6j l39| US] , but the application of the model to experimental data has been hampered by an inability 
to determine the age-dependent model parameters. The usual approach for parameter estimation is 
to solve numerically an inverse problem (see [D El El UHl HH [22l EH EHl SOI SD 1121 H31 [37] on this 
question for structured population models), but this requires extensive input data and is specific to 
a given situation. A much more convenient approach is to assume that the distributed parameters 
lie in a class of functions with only a few constants (power laws or other forms |23j ) and obtain the 
parameters from the fit of the functions to experimental data. 

One function that often provides a better fit to IMT data than others, such as log-normal, inverse 
normal or gamma functions, is an exponentially modified Gaussian (EMG) [251 144] . Under conditions 
in which the IMT distribution can be explained by an EMG model, we submit that the age-dependent 
division rate can be identified as an error function. Starting from this observation, we present a simple 
method to recover the parameters of this error function from parameters fitted to the experimental 
IMT data. Once reliably parameterized, the age-structured model can be used to make predictions 
about cell age-dependent effects of perturbations, for example, whether cells arrest during their cell 
cycle in response to treatment with antiproliferative compounds. 



2 Intermitotic time and age-structured model 

As with human populations, we can associate an age to each individual cell in a cell population. 
We define the age of a cell as the time elapsed from its last mitosis. Thus, the IMT of a cell is its 
age at division. This definition allows us to interpret the IMT distributions in terms of a dynamic 
age-structured population model. 

Experimental IMT distributions can be seen as histograms which represent, for a given population, 
the density of cells with a certain age of division (see Figure Q] for an example taken from [44] ) . The 
age of division is distributed into N bins of width Aa. For all i between 1 and N, the height Hi of 
the i th bar represents the density of cells with an age of division in the window [iAa, (i + l)Aa]. The 
histogram (i?j)i<j<Ar is normalized to represent a density 



We briefly explain here the method used in [44] to build IMT histograms (we refer to this paper 
for more details). The data are obtained using extended temporally resolved automated microscopy 
(ETRAM) in which cell nuclei are fluorescently labeled, imaged by automated time lapse fluorescence 
microscopy, and tracked as individual cells from the resultant image stacks. Cells are subjected to 
various microenvironmental conditions (such as the addition of a drug) at a specific time during image 
acquisition (set to time to) and the effect of the perturbation from that point in time is followed across 
the entire population and individual cells within it. The duration of observation (T) is chosen large 
enough to observe that almost all cycling cells divide before this final time. The data are organized 
into the bins to obtain a histogram, which is then normalized to represent a density of cells. 
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Figure 1: Plot of an histogram (i^j)i<i<jv representing the IMT distribution of a population of PC-9 
lung cancer cells. The data is from [33] and the bin width is Aa = 



The next step is to describe the IMT distributions in terms of an age-structured model and use this 
interpretation to parameterize the renewal equation, also known as the McKendrick-Von Foerster's 
model, which is widely used to model the cell cycle. The partial differential equation for this model 
provides the evolution of the density p(t, a) of cells with age a at time t 



d d 

— p(t, a) + g^p(t, a) + /3(a)p(t, a) 



0, 



t > 0, a > 0, 



p(t,0) 
[ p(0,o) 



2 / /3(a)p(t,a)da, 
/o 



(2) 



Po(a). 



In this model, the cells age one-to-one with time at speed % = 1, and divide with rate /3(a) > 0. 
When a cell divides at mitosis, it produces two daughters with age a = 0, which is taken into account 
by the boundary condition at a = 0. The number of cells at time t with age between a\ and 02 is 
fai °) an< ^ ^ ne total number of cells at time t is / °° p(t, a) da 

The division rate /3 has a probabilistic interpretation: the probability that a cell did not divide by 
age a is given by 

P(a) = e -Jo a 

In this basic form of the model all cells are considered to be proliferating, and all must divide at some 
time (models including cell cycle arrest and death are described in subsequent sections). Thus, the 
division rate has to satisfy 



lim 

a— >+oo 



/3(a') da' = +00. 



(3) 



Equation (|2"D is a transport equation which satisfies a maximum principle, namely if the initial 
distribution po(a) is nonnegative (positive) then the distribution p(t, a) remains nonnegative (positive) 
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for all time t > 0. The solutions of (2) have a remarkable behavior as time evolves, in that the age 
structure equilibrates no matter what the initial age structure of the population may be. This effect 
is known as asynchronous exponential growth, and its interpretation is that the population disperses 
over time to a limiting asymptotic equilibrium age structure where the fraction of the population in 
any age range [01,02] satisfies 

f£p{t,a)da 
t^oo J °° p(t, a) da 



= a constant independent of the initial age structure [56]. Moreover, the solutions to this equation, as 
t — > 00, are known to behave like a separated variables solution, that is, p(t, a) = a function of time 
only x a function of age only. More precisely, consider the eigenvalue problem 

Xp(a) + d a p(a) + P(a)p(a) = 0, 
< p(0) = 2j™l3(a)p(a)da, (4) 
, P(0 > 0, Jp(a)da = 1, 
This problem has a unique solution given by 

p(a)=p(0)e- f ^ (a ' )+X)da ' 
where A > is the unique value such that 

poo 

1 = 2/ /3( a )e-/oW)+AMa' da (5) 

J 

and 



£(0)= QfV/oWHA)*.'^ 1 



Then we can prove that, for large times, 

p(t,a) ~ constp(a)e xt 

(see Appendix lAl for more details and references). If the population of cells proliferates over a suffi- 
ciently long time, we can assume that this asymptotic behavior is reached and use it to investigate 
the IMT distributions. An experimental observation that the total population (independent of age 
structure) is growing exponentially is an indicator that the population has effectively reached the 
equilibrium age distribution, which can be checked using ETRAM for other time series data collection 
(see Figure [2]). 

Now we give a continuous expression of the IMT distribution in terms of the age-structured model. 
The age distribution of the cells relative to time to (time of perturbation) is given by a truncation of 
the equilibrium age distribution p 



Po(a) 



p e - Jo (/V )+ A ) da ' if0<a<t , 
if a > to- 



where p is a scaling constant. We then follow this age distribution along time and obtain for t > 

_ u , ( -tf(P(a>)+\)da> e \t if t <a<t + t , 

P(t,a) = \ n ~ .- 

U otherwise. 
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Figure 2: The evolution of the total population N(t) is plotted on log-scale (blue) (the data is 
from |44j ) . The graph is well-fitted with a line of slope A = 0.022 (correlation coefficient: R 2 = 0.9967). 
This gives evidences that the total population is growing exponentially fast with exponential constant 
A. Thus, N(t) = exp Ai iV(0), and if t* is the population doubling time, then N(t*) = 2N(0), and 
i* = ln2/A. 

According to the age-structured model, and because of the normalization ([1]), the IMT distribution 
satisfies 

J T (o) := C7f 1 / /3(a)p(t,a)dt, C T := / / (3(a)p{t,a) dtda, (6) 
Jo Jo Jo 

by definition of the division rate /3. The fact that no cell can divide in a time less than to means 

Va<t , /3(a) = 0. (7) 
Under this condition, for T large, the function It is close to the function 

loo (a) :=/3(a)e-tf /3 (°'> da '. (8) 

This convergence is made mathematically precise by the following claim (the proof is given in Ap- 
pendix[B]for a more general case in which it is not assumed that p(t, a) is close to the equilibrium age 
distribution): 

Claim 1. Under Assumptions (|3|) and (J7J), we have the convergence 

|/ T (a) - Ioo(a)\ da > 0. 

T— s-oo 

Remark 2. We can easily prove that, under the additional assumption 

lim /3(a) e -Jo a ^ a ')rf»' =0 , (9) 

a— >-+oo 
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we also have 



sup \It(o) - loo (a) | > 0. 

a>0 T^+oo 



Condition ([9]) is satisfied from Assumption ([3]), if for example, f3 bounded or monotonic. 

Because in the experimental protocol T is chosen large enough to observe no dividing cells at the 
end, we can approximate It by which has a simple expression in terms of the division rate ([8]). 



3 From the intermitotic time to the division rate 

The expression Iqq can be inverted to recover the division rate from the IJVIT distribution (see |3Q[ fl3] ) 

/3(a) = rJrrkr-r ^ 

J a Ioo(a')da' 

We start from the fitting of the experimental IMT distributions, which are observed to be positively 
skewed. Prom a fitting procedure for 1^, we use Equation fjlQf) to recover the division rate of the 
renewal equation, which is very important for applications of age-structured models (see Section [5]). 

First consider as in |30] that the IMT distribution is a shifted gamma function (see also |12j and 
references therein). Setting 

f if < a < m, 

Too (dm, a) = < a — m a-m (11) 
k — e " if a > m, 

we can solve explicitly Equation (|10p and we find 

if < a < m, 

«H a ~ m if a > m . (12) 
o"(ct + a — m) 

So by fitting an experimental IMT distribution with a shifted gamma function, we obtain two param- 
eters m and a which allows reconstruction of the division rate of the renewal equation. 

To have a smoother transition at the minimum age of division to, one can consider a second shifted 
gamma function 

if < a < m, 

Ioafalm,^) = { (a-m) 2 _«^m (13) 

-e ct if a > to. 



Then the corresponding j3 is 



/3(a) 



2a 3 



f if < a < to, 



a 2a 2 + 2a (a — m) + (a — to) 2 
The different functions (fTTj) to (|14p are plotted in Figure [3] 



1 (a-m) 2 „ (14) 

if a > to. y 1 
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Figure 3: The two different gamma functions and their corresponding division rate /3(a) are plotted 
for coefficients m = 17 and a = 2. 

It has been observed that an exponentially modified Gaussian (EMG) is often a better model for 
IMT distributions than the gamma function |25[ EH] . An EMG is defined as the convolution of a 
Gaussian with a decreasing exponential, but after solving it can be written with three parameters as 

iMfkw) = !3 Erfc(^)e-M^-^). (15) 

where the (complementary) error function is defined by 

2 f z 2 
Erfc(z) = 1 = / e~* dt. 

Jo 

Replacing 1^ by an EMG in Equation (|10p . we cannot compute explicitly the expression for f3. But by 
numerical comparison, we obtain a division rate /3 that is essentially indistinguishable from an error 
function (see Figured]). 

Instead of fitting the IMT distribution by an EMG and then fitting the corresponding (3 by an error 
function, we may directly assume that (3 is an error function 

(3(a) = (3 Erfc(^^). (16) 

We can then explicitly derive a new fitting formula for 1^ due to Equation (JSJ) 

JooWA^a) =^Erfc(^)e-fo^ErM^L)^ (17) 

where the integral in the exponential can be computed as 

/ Erfcl ] da' = mErfcl — ) =e~ { ^> -(m-a)Erfc( ) + —=e~ { —> . (18) 

7 V a J V a J t/7t \ a J y7r 

Using this formula for the IMT instead of an EMG formula, the fitting parameters provide immediately 
the division rate f3. 
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Figure 4: Left: The division rate /3(a) is obtained numerically from (|10p by using an EMG with 
coefficients m = 22, a = 2 and ft$ = 0.2 to fit 1^ (|15|) . Right: The formula for ft, as an error function 
with these same three parameters, is numerically indistinguishable from the numerically obtained ft 

(R 2 = i)- 



When we know the Malthusian growth parameter (population growth rate A) from experimental 
data (see Figure [2]), it is possible to recover a division rate ft such that relation <|5]) is satisfied. This is 
an important step in order to parameterize and apply an age-structured model to experimental data 
(see Section [5] below) . We notice the relationship of ft in terms of Ioo(a) is 

poo poo 

2 I 00 (a)e- Xa da = 2 ft(a)e~ K ^ a>) da ' e~ Xa da = 1. 
Jo Jo 

To take advantage of this relationship between the IMT data, ft (a), and the malthusian parameter A 
at the equilibrium distribution, we define a new histogram (i?i)i<i<jv by 

Vi, Hi := (19) 

where a-i := (i + ^)Aa is the mean age of the i th bar. Thus defined, the new histogram incorporates 
information about A and satisfies the relation 

oo 

AaJ^fli = l. (20) 
i=i 

We then fit this new histogram (Hi) with the model 

I 00 (a):=2I 00 (a)e- Xa (21) 



instead of fitting (Hi) with Too (a)- Because of the normalization (|2Qj) . we expect that the fitting provides 

/'OC 




parameters such that Ioo( a ) da ~ 1, and this relation can be checked numerically a posteriori. 
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Method and example. We divide the process in three steps and illustrate it by an example. Each 
required fitting step can be performed using the freely available Ezyfit Matlab toolbox 
[< |http : //www . f ast . u-psud . f r/ ezyf it /| >] . 



Step 1: Determine equilibrium IMT distribution 

a) Obtain a histogram (Hi) for the experimental IMT distribution of control (untreated) cells. 
Here we use an IMT distribution obtained in \44^ for PC-9 lung cancer cells using ETRAM 
(see Figure[J\). 

b) Plot the time evolution of the total population on a log e -scale, verify a linear fit and obtain 
the slope as the experimental value for the Malthusian parameter A of the cell population 
(see Figure^ for an example). 

c) Construct the new histogram (Hi) from (Hi) and A by using the definition (|19|) . 

Step 2: Obtain parameters from model fit to IMT distribution 

a) Choose a form for 1^ as a gamma function (jllf) or (I13p . or as the new EMG form (|17p . 

b) Fit the histogram (Hi) with the corresponding form 1^ from definition (|2ip . For PC-9 can- 
cer cells we choose the form (|17l) . because in \44^ it was observed that the IMT distribution 
appeared to be an EMG (see Figure^ for the example). 

c) Verify that the correlation coefficient R 2 of the IMT data and the chosen form Ioo(a) is 
close to 1. 

Step 3: Parameterize age-structured model 

If numerical integral J* °° Ioc(a) da is close to 1, the fitting parameters provide a good approxi- 
mation of the division rate /3, which can then be propagated through the population, based on 
the choice of 1^ . In the example, /3 is given by Equation ([To]) with the numerical paramaters of 
Figure 



4 Model with a death rate 



If it is known that the proliferating cells have a mortality rate (i(a) > 0, then this mortality rate can 
be introduced in the renewal equation: 



-^p(t,a) + —p(t,a) + P(a)p(t,a)+Lt(a) P (t,a) = 0, t>0, a > 0, 



p(t,0) = 2 I3(a)p(t,a)da, 
Jo 

k P(0,o) =p (a). 
We have the same result as before about the long-term behavior 

p(t,a) ~ constp(a)e xt 



(22) 
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Figure 5: Fitting of experimental data in |44j with the model (|17p . The fitting parameters are 
A) = 0.14204, m = 24.456 and a = 3.3451. The correlation coefficient is R 2 = 0.95363 and the 
integral of Ioc(a\(3o, m, a) is f °° Ioo(a) da 1.0983. The formula for the age dependent division rate is 



where the eigenfunction p(a) is now 

p(a)=p(0)e-/oW)H-Ma')+A)^ 

and the eigenvalue A satisfies 

/'OO 

/5(a)e-/o a (/3K)+MK)+A)da' da = 1 



(23) 



Assuming that the stable distribution is reached, we define as in Section [2] the IMT function 



I T (a) := C^ 1 / (3(a)p(t,a) dt 



with 



p{t,a) = pe-^^ + ^ + ^ da 'e x H t < a < t+t0 and C T = / P(a)p(t,a) dtda. 

Jo Jo 



Adapting the proof of Theorem [3] in Appendix |Bl we obtain that It(o>) converges as T — > oo to 

Joo(a) :=C- 1 p(a)e-^(P( a ')+^'))<^' (24) 

where 

poo 

(25) 



Cna '■- 



PWe-KWWV** da . 



For this model there does not exist an inversion formula as Equation (|10p to recover (5 analytically. So 
we consider that /3 is a rational fraction given by Equation (|12|) or (|14|) . or an error function given by 
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Equation (|16p . and we derive the corresponding model 1^ from Equation (|24p . To use this model to fit 
IMT distributions, one has to prescribe the shape of the death rate. If we do not have any information 
about its age dependency, we assume that /i(a) = /i is age-independant, and thus an additional fitting 
parameter (a fixed rate based on measured data can also be used if available). 

The constant Coo was equal to 1 in model ([2]), but if there is a positive death rate, we have that 
Cqo < 1- This constant is a function of the fitting parameters (see Equation (j25fl ) . but we do not have 
an analytic expression of this function in general. A solution to this problem is to use the Malthusian 
parameter A, which incorporates information about \x and /3(a), to define a form ioo(a) as in Section [3] 
Indeed, because of relation (|23p . we obtain 

/00(a) := 2/3(a)e-^ a (/ 3 K)+M)^' e -Aa 

which does not involve Cqo- Then, the division rate /3(a) is obtained from loo (a) as before using the 
modified histogram (Hi) defined from (Hi) by (fl9l) . 

Example. We fit the same distribution as in Section still considering that f3 is an error function. 
With a constant death rate fx, we obtain the four parameters model 

UalA),™,^) = 2p Erfc (^) e~ ^Erfc^)^)^ (26) 

where the integral J Q a /3q Erfc( m ~ a ) da' is given by Equation (|18p . The fitting provides new param- 
eters for f3 and a positive death rate \i (see Figured]). The correlation R 2 is slightly better than in 
Figure [5] and the integral Jq 00 ioo (a) c^a is closer to 1. So we can consider that mortality has to be 
considered for this cell line. 



5 Application 

Once accurately parameterized, model ^) or (|22|) can be used for further investigations into the 
dynamics of age-structured cell populations. We give here one application for in vitro PC-9 cancer 
cell lines undergoing treatment by the drug erlotinib. It has been recently shown in [33] that the 
main effect of erlotinib on cancer cells is to induce entry into quiescence. In [33] a system of ordinary 
differential equations (without age structure) is used to model these experiments. We hypothesize 
here that erlotinib induced quiescence is linked to the age of the cells involved. Many age-structured 
models with quiescence can be found in the literature (see [3 QUI CEH ETJ [28l [29] for examples). Here 
we present a simple age-structured model with quiescence which allows us to explain observed delays 
in response to erlotinib. We start from the observation that there is no effect of the treatment on total 
population growth during the first twenty hours (see Figure [7]). This time corresponds almost exactly 
to the minimal age of division observed for PC-9 cells. It suggests that erlotinib acts only during a 
specific phase of the cell cycle, which based on its biological activity would be expected to be in Gl. 
The model we present is based on this idea and considers a fractional rate / of cells that become 
quiescent in an age-dependent manner, where the fraction is assumed to reflect the dose of erlotinib 
used for the treatment. We estimate the value of / from data available in [33], and then compare 
numerical simulations of the model to the experimental data. Many numerical methods are available 
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Figure 6: Fitting of experimental data with the model (|26p . The fitting parameters are /3q = 0.17879, 
m = 25.007, a = 3.6141 and // = 0.00333. The correlation coefficient is R 2 = 0.95611 and the 
integral of I OQ {a\j3Q, m, a) is f °° ioo(a) da « 1.0132. The formula for the age dependent division rate is 
P{a) = faErf 



for structured population equations (see for instance [2J El El El |2U EI] ) ■ Here we use a scheme based 
on the method of characteristics, as in [HEIET], for its anti-dissipative properties. 



We start from the renewal equation (|22p and we introduce the quantity Q(t) of quiescent cells at 
time t. This quantity evolves according to an ordinary differential equation which is coupled to the 
renewal equation through the system 

—p(t,a) + —p(t,a)+/3(a)p(t,a) + fj,p(t,a) = 0, t > 0, a > 0, 



/•oo 

p(t,0) = 2(1-/) / p(a)p(t,a)da, 
Jo 



p(0,a) =p (a), 
d 



(27) 



it Q(t) = 2fJ P(a)p(t,a)da-vQ(t), 
{ Q(0) = Q . 

The coefficient v is the death rate of quiescent cells and / € [0, 1] represents the fraction of proliferating 
cells which become quiescent. The erlotinib treatment starts at time t = 0, so at this time there are no 
quiescent cells (Qo = 0). The age distribution of the proliferating cells is assumed to be at equilibrium 
as in Equation (f22j) (i.e. po(a) = const p(a)). The experimental values of the total population N(t) 
along time are ploted after normalization by the initial value N(0) on a log-scale (see Figure [7]). 
Because of this normalization, we consider an initial distribution such that J °° po(a) da = 1 which 
leads to po(°0 = p{o) because of the definition of p. 
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The coefficients /3 and /i are parameterized from IMT distributions for cells without treatment using 
the model (|26p . To estimate the value of the fraction / for different doses of treatment, we use data 
available in [44] . The fraction F of quiescent cells is estimated by examination of whether cells treated 
with erlotinib divide or not before the end of the experiment. We want to use this experimental 
fraction F to estimate the coefficient / of the model. For the sake of simplicity, assume that the death 
rates \x and v are in ([27]) . In this case, at the end of the labeling period to, the quantity of labeled 
quiescent cells corresponds to Q(to) and the quantity of proliferating cells corresponds to J *° p(t, 0) dt. 
Thus, the fraction F of quiescent cells at to is 



Q(t ) 
Q(t ) + J to p(t,0)dt 



F= - ; "T,. ~ - ■ (28) 



Now we compute this quantity from model (|27p . keeping in mind that we have assumed no mortality. 
We have 







and 



Finally we obtain 



dt 

/•to roo 

= 2f / p(a)p(t,a)dadt 
Jo Jo 

t rto roo 

p(t,0)dt = 2(1 - /) / / p(a)p(t,a)dadt. 
o Jo Jo 



Q(to) V 



Q(to) + Jo°P(t,0)dt 2/ + 2(1-/) 

So when there is no death rate, the experimental fraction F corresponds exactly to the coefficient 
/. In our simulations we consider positive death rates as suggested in Section [31 Since we cannot 
directly measure the coefficient v, we assume that its value is close to the value of \x. Because the 
numerical values of v and n are similar and very small (/x <C 1), we can consider that / is still well- 
approximated by F and we use the fractions found in [33] to parameterize /. We can see in Figure [7] 
that the numerical simulations are very similar to the experimental curves. In particular, the delay of 
twenty hours before the effect of the treatment on the growth of the population is apparent, and this 
twenty hour period pulses two more times as the population approaches a new equilibrium distribution 
during the total time of the experiment. We have thus developed an age-structured model that can 
explain the dynamic effects of erlotinib on PC-9 cells, which are intrinsically dependent on the age of 
proliferating cells. 



Conclusion 

Linking experimental observation of cell behavior between the single-cell and population scales has 
recently been described using newly developed mathematical models [33] • However, this approach 
does not take into account the possible cell age-dependent effects of a perturbation on cell behavior, 
such as would be expected if the effects occur in cells at a specific position in the cell cycle. Since 
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Figure 7: Left: experimental data without erlotinib, with a 50nM dose and with a 5000nM dose. The 
total quantity N(t) is plotted on a log-scale. Right: numerical simulation of model ([271) for / = 0, 
/ = 0.6 and / = 0.84. The curves represent the evolution of \n(P(t) + Q(t)) where P(t) := J" °° p(t, a) da 
is the total quantity of proliferating cells at time t. They are obtained by solving numerically Equa- 
tion (|27p with the parameters of Figure [6] and with v = 0.004. 

these studies were performed with asynchronously dividing cell populations, it is evident that the 
mathematical models of these experiments should have age structure as a primary feature. In fact, 
to fit the data, the authors had to use an artificial time offset to account for the age-structured 
effects. Here we provide a formalized approach to accurately account for cell age-dependent effects 
on cellular behavior. A major difficulty in the parameterization of age-structured models is the 
determination of the age dependent division rate. Our study provides a method for the quantitative 
recovery of this rate by fitting experimental IMT distributions to special forms, such as gamma 
functions or exponential modified Gaussians. This model, once successfully parameterized, is very 
useful for simulating and analyzing age-dependent phenomena in cell population dynamics. We have 
presented one such application for the in vitro treatment of cancer cells by erlotinib. This example 
shows the utility of age-structured population models in explaining the connection of drug therapy 
to phenomena such as cell cycle phase entry into quiescence. The method we have presented can be 
implemented readily to many issues in cell population behavior when there is experimental data based 
on cell age. 
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A Convergence to the equilibrium 

The long-time behavior p(t,a) ~ constp(a)e xt can be proven by using semi-groups methods (see [341 
145] ) or General Relative Entropy (see [351 [36]) techniques. We detail here a result provided by General 
Relative Entropy. 

Consider <j) the unique solution to the adjoint eigenvalue problem 

J \<j>(a) " da^a) + /3(a)0(a) = 20(O)/3(a), 
\ 0(0 > 0, / p(a)(p(a)da = 1. 
Then the General Relative Entropy method (see [35[ [36] ) allows us to prove that 



where 



\p(t, a)e xt - p p(a)\(p(a) da >0 

t— >oo 



oo 

Po = I <f>(a)p (a) da. 
o 



B Convergence of the IMT distribution 

Theorem 3. Suppose that Assumptions ([3]) and ([7]) are satisfied and consider an initial distribution 
po such that po(a) = for all a > to- Then we have the convergence 



/•oo 

/ \I T (a)-Ioo(a)\da >0 

Jo T ^°° 



where It and 1^ are defined in ([6]) and ([8]). 

Proof. First we extend the initial distribution and the division rate /3(a) to negative ages by setting 
Po(a) = (3(a) = for a < 0. Since the daughters of the labeled cells are not tracked, the age distribution 
p(t, a) of labeled cells satisfies a transport equation without boundary condition. Thus it writes, using 
the characteristic method, 

Vt,a> 0, p(t,a)=po{a-t)e-ti p{ - a - s)ds . 
Introducing this expression in the definition of It{o) we obtain, with changes of variables, 

r-T 

I T {a) = C^ 1 / (3{a)p {a - t)e~ ^ 
Jo 

= Cy 1 f f3(a)p (a-t)e-^-^ {a ' )da 'dt 
Jo 

= Cy 1 ^ ^{a)p {u)e"^^ da 'du 



a-T 

a 



= C^ 1 ^{a)e-^^ a ') da ' / po{u)e^^ da 'du. 

Ja-T 

Since the support of po is included in [0,io] and /3(a) = for a G [0, io] due to Assumption ([3]), we 
have for all u £ R 

p (u)e^^ da ' =p (u) 
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so It (a) writes 

I T (a) = Cy 1 /3(o)e-io a/3(a ' )da ' / p (u)du. 



a-T 



Define the primitive 

P(a) := / po(u)du. 



This function is nondecreasing and satifies P(a) = for a < and P(a) = P(to) = Jq°Po f° r a > 
so we have 



I p (u)du = P(a) -P(a-T) = < 

J a-T 



' P(a) if < a < t , 

P(to) tft <a<T, 

P{t )-P(a-T) HT<a<T + t , 

if a > T + t . 



Because /3(a) = for a < to, we obtain by integration of It{o) on R + 

pT fT+t 

l = C^P(t ) / p{a)e-Kft a ' )da ' da + C? 1 f3(a)e~ So ^( a ') da ' (P(t ) - P(a -T))da 

Jo Jt 

which gives 

pT pT+t 

C T -P(t ) / f3{a)e-^^ da ' da= / /3(a)e~ •fo 0( a ') da> (P(t ) — P(a — T))da 
Jo Jt 

T+t 



pi +to 

<P{t ) J (3{a)e-^^ a ' )da ' da 



T^+oo 



-> 0. 



Since 



we conclude that 



poo 

/ j3{a)e-^^ a ' )da ' da = 1, 
Jo 

poo 

C T ► P(t ) = / Po(a) da. 

T^+oo J 



Then we have 

coo r T 



poo pi 

/ |/t(o) - Joo(o)| da < \C^P(t ) - 1| / ^(^e-io"/ 3 ^')^' da 
Jo Jo 

pT+to 

+ y \C^[P(t ) -P(a-T)} - l\l3(a)e-foP( a ') da ' da 



and it ends the proof of Theorem [3j 



+ / f3{a)e-^ ^ da ' da 

< \C^P(t ) - l| / p(a)e-ftK a ' )da ' da + K / p{a)e~ So P( a ') daf da 
-> 



o 
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